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Abstract 

Electrical stimulation (ES) devices interact with excitable neural tissue toward eliciting action potentials (AP's) by specific 
current patterns. Low-energy ES prevents tissue damage and loss of specificity. Hence to identify optimal stimulation- 
current waveforms is a relevant problem, whose solution may have significant impact on the related medical (e.g. 
minimized side-effects) and engineering (e.g. maximized battery-life) efficiency. This has typically been addressed by 
simulation (of a given excitable-tissue model) and iterative numerical optimization with hard discontinuous constraints - e.g. 
AP's are all-or-none phenomena. Such approach is computationally expensive, while the solution is uncertain - e.g. may 
converge to local-only energy-minima and be model-specific. We exploit the Least-Action Principle (LAP). First, we derive in 
closed form the general template of the membrane-potential's temporal trajectory, which minimizes the ES energy integral 
over time and over any space-clamp ionic current model. From the given model we then obtain the specific energy-efficient 
current waveform, which is demonstrated to be globally optimal. The solution is model-independent by construction. We 
illustrate the approach by a broad set of example situations with some of the most popular ionic current models from the 
literature. The proposed approach may result in the significant improvement of solution efficiency: cumbersome and 
uncertain iteration is replaced by a single quadrature of a system of ordinary differential equations. The approach is further 
validated by enabling a general comparison to the conventional simulation and optimization results from the literature, 
including one of our own, based on finite-horizon optimal control. Applying the LAP also resulted in a number of general ES 
optimality principles. One such succinct observation is that ES with long pulse durations is much more sensitive to the 
pulse's shape whereas a rectangular pulse is most frequently optimal for short pulse durations. 

Citation: Krouchev Nl, Danner SM, Vinet A, Rattay F, Sawan M (2014) Energy-Optimal Electrical-Stimulation Pulses Shaped by the Least-Action Principle. PLoS 
ONE 9(3): e90480. doi:10.1371/journal.pone.0090480 

Editor: Dante R. Chialvo, National Research & Technology Council, Argentina 
Received November 11, 2013; Accepted January 30, 2014; Published March 13, 2014 

Copyright: © 2014 Krouchev et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: The Fonds de recherche du Quebec - Nature et technologies and the Natural Sciences and Engineering Research Council of Canada provided funding 
for this work. S.D. was supported by the Vienna Science and Technology Fund (WWTF), Proj.Nr. LS1 1-057. The funders had no role in study design, data collection 
and analysis, decision to publish, or preparation of the manuscript. 

Competing Interests: The authors have declared that no competing interests exist. 

* E-mail: Nedialko. Krouchev@polymtl.ca 



Introduction 

Electrical stimulation (ES) today is an industry worth in excess 
of 3 G$. ES devices interact with living tissues toward repairing, 
restoring or substituting normal sensory or motor function [1]. 
The rehabilitation-engineering applications scope is constantly 
growing: from intelligent limb prosthetics and deep-brain stimu- 
lation (DBS) to bi-directional brain-machine interfaces (BMI), 
which are no longer just about recording brain activity, but have 
also recently used ES toward closed-loop systems, [2-5]. 

Application-specific current patterns need to be injected toward 
reliably eliciting action potentials (AP's) in target excitable neural 
tissue. To prevent tissue damage or loss of functional specificity, 
the employed current waveforms need to be efficient. This may 
significandy impact the biomedical effects and engineering 
feasibility. Hence, an optimization problem of high relevance to 
the design of viable ES devices is to minimize the energy required 
by the stimulation waveforms, while maintaining their capacity for 
AP triggering toward achieving the targeted functional effects. 

A number of recent studies of ES optimality are based on 
extensive model simulation and related numerical methods 



through the wider spread of high-performance computing, e.g. 
[6-9]. The model dynamics to iterate can be arbitrarily complex 
and nonlinear. This implies lengthy numerically-intensive compu- 
tation, irregular convergence and constraints that may be difficult 
to enforce - e.g. that an AP is an all-or-none phenomenon. Thus, 
any function of membrane voltage will suffer dramatic disconti- 
nuities at parameter-space manifold boundaries where intermit- 
tent AP's are likely to be elicited. 

Hence, such an iterative approach is not only computationally 
expensive, but its solution quality is highly uncertain and model- 
specific. The long-lasting iteration may converge to shallow local 
energy-minima. Such numerical misdemeanor of the approach is 
well known to its frequent users. 

In this work we follow the ES pioneers - we use physical 
reasoning and related mathematics toward a more theoretical 
treatment of the subject. 

Below we summarize very briefly our historical premises. ES' 
theoretical cornerstones were laid a century ago by experimental- 
ly-driven assumptions and models, [10-12]. Various constant ES 
current levels and durations were tried systematically. E.g. Louis 
and Marcelle Lapicque spent many years performing such lab 
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experiments with multiple physiological preparations [13,14]. This 
classical work led to concepts like strength-duration curve (SD), i.e. the 
function of threshold (but still AP-evoking) ES current strength on 
duration. The first mathematical fit to this empirical results is 
usually attributed to Weiss, [10,15] 

lTHR(T) = b(\+c/T) (1) 

where T is the stimulus duration, b is called the rheobase (or 
rheobasic current level) and c is the chronaxie. 

The most expedite way of introducing the rheobase and chronaxie 
would be to point to eqn. (1) and notice that: 

lim I THR (T) = b (2) 

and 

I T HR(c) = 2b (3) 

i.e. the rheobase is the threshold current strength with very long 
duration, and chronaxie is the duration with twice the rheobasic 
current level. In the pioneering studies electrical stimulation was 
done with extracellular electrodes. 

Eqn. (1) is the most simplistic of the 2 'simple' mathematical 
descriptors of the dependence of current strength on duration, and 
leads to Weiss' linear charge-transfer progression with T, 
Q(T) = T x IxHR = b x (T + c). Both Lapicque's own writings - 
[11-13], and more recent work are at odds with the linear-charge 
approximation. Already in 1907 Lapicque was using a linear first- 
order approximation of the cell membrane, modeled as a single- 
RC equivalent circuit with fixed threshold: 

b be- Tit 

Ithr{T) = ^-^777 = b + ^^777 (4) 

with time constant x = C/g; C and g=l/R are the membrane 
capacity and conductance respectively. 

The second form of eqn. (4) is easily obtained by subtracting/ 
adding the term be~ T l x . From it, when T»T (and hence 
e-^^l): 

lTHR(T)*b(l+T/T) 

which accounts for the hyperbolic shape of the classic Lapicque 
SD curve. 

Originally, eqn. (4) described the SD relationship for extra- 
cellular applied current. However, the single-RC equivalent circuit 
with fixed threshold, where I is the electrode current flowing 
across the cell membrane 

Cv+v/R = I (5) 

can be used with either extra- or intra-cellular stimulation. 
v = ( V — V rest ) is the reduced membrane voltage with V rest the resting 
value of V . From eqns. (4) and (5), one may also see that 
b = g(VrHR— V rest ), where Vjhr is the attained membrane 
voltage at the end of the stimulation (at time 7"). 

Notice that the chronaxie c is not explicitly present in eqn. (4). 
Notice also that - with very short duration T«x, by the Taylor 
series decomposition of the exponent (around 7" = 0), one may 
have either Ithr(T) Kbz/T or I T HR(T) = b[\ + x/T\. Note that 
these two different simplifications (and esp. the latter) are 



'historical' and depend on which of the two right-hand sides 
(RHS') of eqn. (4) is used. In the second case only the denominator 
is developed to first order, while the numerator is truncated at 
zero-order. The second approximation throws a bridge to Weiss' 
empirical formula of eqn. (1). I.e. the latter is a simplification of a 
simplification (i.e. of the lst-order linear membrane model), 
capturing best the cases of shortest duration. On the other hand, 
IrHR(T)xbx/T leads to a constant-charge approximation. 
Interestingly, the latter may fit well also more complex models 
of the excitable membrane, which take into account ion-channel 
gating mechanisms, as well as intracellular current flow, which 
may be the main contributors for deviations from both simple 
formulas. These 'subtleties' are all clearly described in Lapicque's 
work, but less clearly by one of the most recent accounts in [16]. 

Before we continue, it is in order to examine the practical value 
of numerical optimization to identify energy-efficient waveforms. 
It is limited for the following reasons. First, it is subject to the 
rigorous constraints of quantitative equivalence between the model 
used and the real preparation to which the results should apply. A 
noteworthy example is provided by the very practice of numerical 
simulations: often a minute change in parameters precludes the 
use of a. just computed waveform, which is no longer able to elicit an 
AP in the targeted excitable model. Alas, the same or similar 
applies hundredfold to the real ES practice. 

Second, in the search for minimum-energy waveforms, using 
numerical mathematical programming algorithms, there is no 
guarantee about obtaining a globally optimal solution. 

Finally, such an approach sheds very little light with respect to 
the major forces that are at play, and the key factors which 
determine excitability, such as - for example, the threshold value of 
membrane potential, whose crossing triggers an AP. 

However, the problem at hand is also reminiscent of the search 
for energy-efficiency in many other physical domains - e.g. 
ecological car driving. For centuries, physics has tackled similar 
problems through an approach known as the Least-Action Principle 
(LAP) [17]. 

Thus, we first used simple models to derive key analytical 
results. We then identified generally applicable optimality princi- 
ples. Finally, we demonstrate how these principles apply also to far 
more complex and realistic models and their simulations. 

The modeling and algorithmic part of this work is laid out in the 
next section. First, we introduce a simple and general model 
template. Next we present four most popular specific ionic-current 
models. Each of these can be plugged in the template to describe an 
ES target in a single spatial location in excitable-tissue (or 
alternatively - a space-clamped neural process). 

We then examine the conditions for the existence of a finite 
membrane-voltage threshold for AP initiation. The introduced 
ionic-current model properties are analyzed to gain important 
insights into the solution of the main problem at hand. 

Two very different ways to identify energy-efficient waveforms 
are presented in the last two subsections of the Methods. The first 
relies on a standard numerical optimal-control (OG) approach. The 
second outlines the LAP in its ES form, which is used to derive a 
general analytic solution for the energy-optimal trajectories in time 
of the membrane-potential and stimulation-current. 

The Results section presents the model-specific results, applying 
OC or the LAP. We perform a detailed optimality analysis for 
both the simple and more realistic models. Comparisons between 
the two types of approaches, and the quality of their solutions, are 
made. 

Commonly used abbreviations are summarized in Table 1 and 
symbols - in Table 2. 
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Table 1. 


Commonly used abbreviations. 




Symbol 


Description 


OD 


zero-dimensional, i.e. single-compartment or space clamp 
models; whose spatial extents are confined to a point 


1D 


cable-like, multi-compartment spatial structure; homo-morphic 
to line 


2D etc. 


two- or more dimensional, refers to the number of states that 
describe the excitable system's dynamics 


AIS 


the axon's initial segment 


AP 


Action potential 


ASA 


Adjoint Sensitivity Analysis 


BCI 


brain-computer interface 


BMI 


brain-machine interface 


BVP 


Boundary-value [ODE solution] problem 


BVDP 


the Bonhoeffer-Van der Pol oscillator-dynamics model; also 
known as the Fitzhugh-Nagumo model 


DBS 


Deep-brain stimulation 


ES 


Electrical stimulation 


FHOC 


Finite-Horizon Optimal-Control 


FP 


Fixed point of system dynamics — > vanishing derivative(s) 


HH or HHM 


Hodgkin and Huxley's [model of excitable membranes] 


IM 


the Izhikevich model 


LM 


the Linear sub-threshold model; also known in computational 
neuroscience as leaky integrate & fire 


MRG 


the Mclntyre, Richardson, and Grill model 


OC 


Optimal-Control 


ODE 


Ordinary Differential equation; see also PDE 


PDE 


Differential equation involving partial derivatives; see also ODE 


LAP 


the Least-Action Principle 


RN 


Ranvier-node 


RHS 


right-hand side 


SD 


strength-duration [curve] 


W.R.T. 


with respect to 


doi:1 0.1 371 /journal.pone.0090480.t001 



Methods 

A General Excitability Model Template 

For the equivalent circuit of Fig. 1, Ig is the stimulation current. 
Ic is the capacitive current, whose direction is as shown on the 
Figure when the excitable-membrane's potential is being depolar- 
ized. The algebraic sum of all the ionic and all axial currents is 
represented by 1?:= I ion + 1 axial, where I ax ial stands for the 
algebraic difference (divergence) of in- and out-going axial 
currents. In the sequel we will use the notation u(t) = Is(t) for 
the stimulation-current waveform. The latter is our system input, 
which will be the leverage to refine in order to achieve desirable 
outcome - reliable triggering of APs in the excitable system. It is 
customary in the control literature to denote such a signal u(t). 

Thus, all the currents are linked by the first Kirchhoff circuit 
law: 

u(t) = I s (t) = I c (t) + h [ V(t),x(t)] = C m V + h( V,x) (6) 

where - in the most general form, depends on membrane 
voltage V(t) and on the state vector of the ionic channels' gate 
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Table 2. Commonly used symbols. 





Symbol 


Description 


C or C„ 


membrane capacity 


At 


the temporal precision of a model's simulation 


g or g,„ 


membrane conductance; see also R ni 


gx 


nominal (max.) conductance for ion X 


GE 


the growing-exponent stimulation pulse 


Is 


stimulation current, see also u{t) 


Ic 


the capacitive current, see also C m 


Ithr(T) 


threshold current for duration T to elicit an AP; see 

Tstim 


Iaxial 


algebraic sum of in and out axial currents 


h,n(V(t)) 


ionic current function of membrane voltage; see V(t) 


h,nfi(V) 


resting-state approximation; see xq 


W,(F) 


asymptotic-state approximation; see x c£) (V) 


\ 


cable spatial constant 


R or R,„ 


membrane resistance; see also g m 


P and P(T) 


power for u(t) as function of duration; see u(t), Tstim 


Q and Q(T) 


charge-transfer 


SQR 


square (rectangular) waveform 


Tcit 


critical duration; see Tstim 


Tstim or Ts or T 


duration of stimulation 


t or x m 


membrane time constant 


Tfa> or x x 


gate time constant for ion X 


u(t) 


stimulation waveform 


u'(t) 


optimal current stimulation waveform 


V 


membrane voltage 


V, or V ral 


resting V 


v=v-r R 


voltage difference w.r.t. rest 


V or dVjdt 


first time-derivative of the membrane voltage 


no 


temporal pattern of V 


V(t) 


optimal V(i) 


Vthr 


AP triggering V threshold 


Vthr$i 


/■esf/ng-state V THR 




the asympfot/c-state Vthr 


Xo=x w {V r ) 


gate resting state for ion X; see V r 


x 0 ,(V)=\im l ^ tx> x(t\V) 


gate asymptotic state for ion X 



doi:1 0.1 371 /journal.pone.0090480.t002 



variables. Unless ambiguous, below we will simplify notation by 
writing Iy.(V). 

C m (typically around 1 \xF /cm 2 , [18]) and V(i) (in mV's) are 
the excitable-membrane's capacitance and potential. Equation (6) 
can be rewritten as: 

C m V = u(t)-h(V) (7) 

Clearly according to eqn. (7), an outgoing total ionic current 
opposes the effects of cathodic stimulation, since not all of u(t) is 
employed toward the main goal of maximizing the V(i) growth, 
which the reader may have also already deduced from the 
equivalent circuit of Fig. 1. Conversely, ingoing total current assists 
the effects of stimulation. Hence, in such a case u(t) may be lower 
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Figure 1. Excitability model template: The equivalent circuit represents the simplified electro-dynamics of an excitable membrane. 

Is is the intra-cellular stimulation current. Iq = CV is the capacitive current. The direction of the latter is for a case of depolarizing the membrane's 
voltage (i.e. the inside of the cell wall becoming more positive). The algebraic sum of all the ionic and all axial currents is represented by 
It. = Iion + 1 axial . where stands for the algebraic difference (divergence) of in- and out-going axial currents. 
doi:1 0.1 371 /journal.pone.0090480.g001 



than when it is estimated assuming the absence of membrane 
conductivity. Let us elucidate right away by providing typical 
examples. 

Specific Single-compartment (Space-clamp) Models 

The models here are zero-dimensional (OD). Their spatial extents 
are confined to a point. This may be contrasted to the multi- 
compartment cable-like models that we will discuss later, and 
whose spatial structure is one-dimensional (ID) - i.e. homo-morphic 
to a line. 

For single-compartment models there are no axial currents. 
Hence, h=IiON- 

Linear Sub-threshold model (LM). 

hoN{V) = g m (V(t)-V r ) (8) 

g m is the excitable-membrane's resting (V= V r ~ — 70 mV) 
conductance - in milli-Siemens per unit membrane surface area 
- e.g. 1 mS/cm 2 . Substituting Iion(V) from eqn. (8) into eqn. (6) 
yields a linear first-order model with %=C m /g m = R m C m the 
familiar expression for the time constant of such a dynamic model. 
This model predicts a reasonable resting x x 1 ins. 

As pointed out in the introduction, this type of model 
was extensively used by the ES pioneers, [12]. They 
were particularly concerned with the derivation of analytic 
expressions for the experimentally observed strength-duration 
(SD) curves. The latter describe the threshold (minimal) 
current strength (Ithr), which if maintained constant (i.e. 
through a rectangular waveform) for a given duration T is 



likely to elicit an AP in excitable-tissue (see the introductory 
section). 

Even if it may account for a significant part of the sub-threshold 
variation of the membrane's potential, the linear model lacks a 
paramount feature - it cannot fire AP's as the latter are due to the 
highly nonlinear properties of the excitable-membrane's conduc- 
tance around and beyond the firing threshold. 

The Hodgkin-Huxley-type model (HHM). Hodgkin and 
Huxley (HH) not only proposed a novel way to model ionic- 
channels but also introduced ionic-channel-specific parameters to 
fit experimental data [19], Since, HH-type models have been 
proposed for many ionic-channels for cardiac to neuroscience 
applications. 

We present one such model from the literature - [20], which has 
been used to fit experimental data from the central nervous system 
and particularly the neocortex. 

Iion( V,x) = g Na m i h{ V - E Nu ) + g K n{ V-E K ) 

+ gleak(V — Eleak) 

See Tables 3 and 4, which define all the model's variables and 
parameter values. We consider specifically the Na^ 6 sodium 
channel subtype, to which the axon initial segment (AIS) owes its 
higher excitability [20,21]. 

The dynamics of a gate-state variable x(t) (where x(i) stands for 
one of m(t),h(t),n(t)) are described by: 
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T x (V)x + x = x m (V) (10) 

Eqns. (6), (9) and (10) define a system of four coupled ODE's - 
with respect to the four dynamic variables [V ,m,h,ri\(i). 

Further simplification may reduce the model complexity, 
maintaining only V(t) as the single dynamic variable. Gate- 
variable states are factored out by introducing appropriate non- 
dynamic functions of the membrane potential. E.g. in eqn. (9), the 
fast m gates may be assumed to reach instantaneously m m (V), 
while the far slower h and n gates remain at their resting values 
(corresponding to a membrane at its resting equilibrium potential 
Vr). 

The Izhikevich model (IM). 

Iion{ V,w) = w - 0.04 V 1 - 5 V - 1 40 (11) 

This model [22] has a second-order nonlinearity, compared to its 
predecessor - the BVDP model [23], which contains a cubic 
nonlinearity. The IM will therefore not auto-limit. As in the 
BVDP, there is a slow second dynamic variable w(t) called the 
'recovery current' and its dynamics is described by: 

w/c = bv — w (12) 

The IM responds to supra-threshold stimulation with a wide 
variety of AP-firing patterns, depending on the particular choices 
of parameters. Interested in the sub-threshold regimen, we have 
chosen the "Spike Latency" set: 6 = 0.2,c = 0.02 [24]. Hence, 



X v = 1 jc is equal to 50 ms. At the time-scale of a single 
stimulation pulse (lasting at most a few milliseconds), W is virtually 
a constant. 

Here, it may be important to remind the reader that the 
state of simplest models like the IM needs to be artificially reset 
after an AP event. However in more complex models (e.g. the 
HHM), channels that are responsible to revert the system to 
its resting potential will have a significant effect on the 
optimal waveform. We will see this in more detail in the results 
section. 

Multi-compartment Models 

To expand the scope of our analysis and the applicability of 
its results, it is essential to also address models of AP initiation 
and propagation along spatial neural structures. A popular 
example is the Mclntyre, Richardson, and Grill model 
(MRG'02). It was originally used to simulate the effects of ES 
in the peripheral nervous system and specifically the myelin- 
ated axons that form nerve bundles [25]. An adapted version 
of the same model was recently used to simulate the effects of 
DBS [7]. 

Myelinated axon has been pinpointed as the most excitable 
tissue with extracellular stimulation [26-28]. Therefore models 
like the MRG'02 are of particular interest. Moreover, this 
model facilitates the illustration of optimality principles as it 
has only one excitable compartment type - the Ranvier-nodes 
(RN). The paranodal and other compartments that form 
the myelinated internodal sections are all modeled as a 
passive double-cable (due to the myelin sheath that 
insulates the extracellular periaxonal space) structure, see 
Fig. 2. 



Table 3. Definition and notation for the key HHM variables. 



Notation 


Variable description and units 


Typical value (*1 


Potentials, in mV: 


V m 


Membrane voltage 


(*3 


V rm 


Membrane resting voltage 


-77 


Ek 


K + Nernst potential 


-90 




Na + Nernst potential 


60.0 


Etaik 


Leak reversal potential 


-70 


Membrane capacitance, in /.iF/cm 2 : 


c 


Membrane capacitance 


1 


Maximum (*2 conductances, in 


mS/cnr: 




SA' 


K + conductance 


150 


gm 


Na + conductance 


300 


gLeak 


Leak conductance 


0.033 


Currents, in fiA/cm 2 : 


Ik 


K + Ionic Current (*4 


gK><nx(V m -E K ) 


Inu 


Na + Ionic Current 


gNaXni s hx(V m -E m ) 


iLwk 


Leak Current 


gL„„k^(Vm-E L „ k ) 


Notes: 



(*1 Typical values are for the Na v \^ model, [20]; see also Table 4. 

(*2 These are dependent on (grow with) temperature, the values listed are for T = 23 C. 

(*3 Membrane voltage is either at its resting value V rcsl ; is depolarized (grows due to stimulation and/or activated sodium Na + ion channels); is repolarized (decays 
back to V rm , due to the potassium A' + ion channels). 

(*4 Ionic currents depend on both the membrane voltage and the dynamic state of the ion channels' gates. See Table 4. 
doi:1 0.1 371 /joumal.pone.0090480.t003 
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Table 4. Gate-state dynamics parameters. 





Notation 


Variable description 


Value 


Temperature dependence: 


gio 


f2io constant (*1 


2.3 


K + : n-gate (*2 


a u 


/?-gate max opening rate 


0.02 


b„ 


//-gate min closing rate 


0.002 


V*M2 


half-min/max in/activation rate voltage 


25 mV 


k„ 


n-gate voltage constant k 


9 


Na+j 6 : ;??-gate (*2 


a m 


m-gate max opening rate 


0.182 


b m 


m-gate min closing rate 


0.124 




half-min/max in/activation rate voltage 


41 mV 


km 


m-gate voltage constant k 


6 


Na+ i : //-gate (*2 




//-gate max opening rate 


0.024 


h 


//-gate min closing rate 


0.0091 


VkM2,a 


half-max activation rate voltage 


48 mV 


^,1/2,6 


half-min inactivation rate voltage 


73 mV 


h 


//-gate voltage constant k 


5 


fc 4 ,» (*3 


asymptotic gate-state voltage constant k c(J 


6.2 


P*,l/2,oo 


50% open gates voltage 


70 mV 



Notes: 

(*1 Temperature dependence is linear and with a slope kr = Q^~ , where 7o = 23"C. 

(*2 For a given gate type y of the K + and /Vfidj 6 ionic channels, the fractions of open and closed gates are given by the general (Boltzmann-Energy like) template 
formulae: 



w=V m -V u 



a y (w) = a y w/(l - e -" /k >') p y (w)= -b y w/(\ -e" ,/t > ) where 
Thus, the corresponding rates of opening doty/dw and closing dft y /dw are sigmoidal functions of w s.t. 

lim a,, = lim B v = 0 lim a. v =a v lim B v = —b v 

The actual position of the inflection point (ir = 0} is determined by the V\j 2 parameter. For the m and n gates, by the I'Hospital-Bernoulli rule, it can be seen that at 
Vm = V\i2< the opening or closing rates attain half of their max or min, respectively. 
(*3 For the inactivating gate /; of the Na^ 6 ionic channel: 

MF) = 1/(1 +e'"'' /lc i^ ) nv, = V m - F WA „ 

doi:1 0.1 371 /journal.pone.0090480.t004 

The RN compartment is a model of the HH-type: 

honi V,X) = gNafm^H V-E Nll )+ g N a,pP 3 (V ~ E Na ) 



+ 9 K n\ V — E K ) + ai eak ( V - E, eak ) 



(13) 



inactivating. In addition, this model has very slow s gates, associated 
to its K + ion channel and very fast m gates. 

Below we call a fixed point (FP) every Vpp value s.t. 
Iion(Vfp) = 0- From eqn. (7) with u = 0, 



= 0 



Here two different Na + ion channel subtypes are modeled 
(please see Table 5 for all the details). The fast subtype (with 
maximum conductance parameter gpfaj) is controlled by the 
opening m and closing h gate states. The persistent subtype (with 
maximum conductance gNa.p) is controlled by the p gates. As its 
name suggests, it has no gate-inactivating states and is non- 



The nonlinear dynamics behavior of the RN compartment 
taken in isolation is quite unlike that of the specific single- 
compartment HHM example we provided above. None of its four 
FPs are stable. Around its unstable 'resting' state [V r = — 80 mV), 
the zero-dimensional RN's of MRG'02 model yield depolarizing 
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Figure 2. The MRG'02 myelinated axon model (See also Table 4) Box: Equivalent circuit for current injection into the center RN 

(#D. 
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ionic current. I.e. not only does I ion not resist moving away from 
the resting state, but it actually contributes to automatic firing, 
with or without any external current! 

The addition of the passive myelinated spatial structures around 
the RN's makes the resting state stable, and the problem at hand 
(of identifying the LAP-optimal ES waveforms) tractable only 
within a spatial structure. However, this also comes with bonuses. 
First, the active-passive association brings a very clear-cut picture 
of the factors at hand that influence AP initiation and propagation. 
Second, the myelinated double-cable has a very low spatial 
constant, which provides for a straightforward extension of the 
single-compartment analysis. 

Namely, consider the second term in the more general 
expression for /j; = I ion + hxial m eqn. (7). Since around the 
resting state I/on is always there as a depolarizing factor, it is I ax ial 
that needs to be closely considered, see Box in Fig. 2. 

The numerical results presented for the MRG'02 in the 
literature [7,8] often target the mid-cable (center) RN in their ES 
simulations. This motivated us to use of the method of mirrors to 
double the model's dimensions at the same computational cost. 
We consider a long axon (with 41 RN's), which has a relatively 
low length constant (J?^/(gaP a ))- See also Tables 4 and 5. For the 
RN's 1= 167.5 jim vs respectively 2129.7 and 443.2 fim, for the 
myelinated and the MYSA (paranode) sections. These are paired 
to significant differences in the passive membrane time constant 
(t ~ Cilg a ). For the RN's t = 0.29 ms vs respectively 20 and 2 ms, 
for the myelinated and paranode sections. The cable end-conditions 
are formed by virtual compartments with membrane at rest 
V r = — 80 mV. This choice is further motivated by the results of 
model simulations - namely the relatively little spread of 
potentials at the end of stimulation lasting up to a few 
milliseconds (see Fig. 3). 



We studied extensively all the published accounts of 
the MRG'02 model and its use for ES modeling [7,8,25]. We 
also carefully compared parameter values (see Tables 5 and 6) 
to the ones in the official NEURON models database 
(senselab.med.yale.edu/ modeldb/ShowModel.asp?mo- 
del = 3810). 

Our model implementation originally for [29,30] was done 
in Matlab (the Mathworks, ver. 7 and above). The code uses 
CVODES (the Lawrence Livermore National Laboratory, 
Release 2.7.0) to reliably and robustly solve the related multi- 
dimensional system of ODEs. The implementation was 
validated through extensive comparisons and personal corre- 
spondence with the authors of the original model - W.M. 
Grill [31] and A.G. Richardson, regarding specifically the 
mismatch between the 2002 publication and its NEURON 
implementation. 

Preliminary Analysis: On the Existence of the AP-firing 
Threshold 

The above ionic-current descriptions differ largely in form and 
complexity. Yet each of them is capable of capturing some of the 
essential dynamics properties of excitable living tissues. 

In order to elicit an AP through electric stimulation, the 
membrane's potential V(t) needs to first be driven (depolarized, 
V>0) to some threshold value Vthr, beyond which assisting ionic 
channels are massively engaged to produce the AP upstroke 
without the need of any further ES intervention. From eqn. (6) in 
order to do so, the stimulation waveform needs to be positive and 
superior to Iz(V,x) at most times - i.e. u(t) needs to overcome the 
opposing currents. 

A Vthr value is hiding inside each of the above nonlinear 
flavors of Ix(V,x). Predictably, it is easiest to find the Vthr value 
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Table 5. MRG'02 double-cable model-axon electrical parameters. 





Notation 


Parameter description 


Value 


Shared parameters: 




Resting potential 


-80 mV 


Pa 


Axoplasmic resistivity 


70 0 


Pp 


Periaxonal resistivity 


70 n 


Nodal compartments: 


c„ 


Membrane capacitance 


2 fiF/cm 


Ek 


K + Nernst potential 


-90 mV 


Enu 


Na + Nernst potential 


50.0 m V 


ELeak 


Leak reversal potential 


-90 mV 




Maximum slow K + conductance with opening .s and no closing gate states 


0.08 S/cm 1 


gm<j 


Maximum fast Na + conductance with opening m and closing h gate states 


3.0 S/cm 2 


gN«.n 


Maximum persistent Na + conductance with opening p and no closing gate states 


0.01 S/cm 1 


gLeak 


Leak conductance 


0.007 S/cm 2 


Internodal compartments: 




Membrane capacitance 


2 fiF/cm 2 


E Ps , 


Passive-compartment Nernst potential 




Passive (leak) membrane conductance by segment type: 


S„ 


MYSA 


0.001 S/cm 2 


Sf 


FLUT 


0.0001 S/cm 2 


gi 


STIN 


0.0001 S/cm 2 


Myelin parameters: 


C>nv 


Capacitance 


0.1 fiF/cm 2 


gmy 


Conductance 


0.001 S/cm 2 


Notes: 

See also Table 6. 

doi:1 0.1 371 /journal.pone.0090480.t005 



associated with the IM. Above we saw that the variable w in the 
IM reacts .slowly to changes in V. Hence, one may approximate it 
by its value at rest: w r = bV r . The resting membrane potential V r is 
then obtained from the condition Iu m fi{V r ) = Q, where the 
subscript 0 indicates that we have assumed w(t) = w r . 

The resting potential V r is one of the zeroes of the 2nd-order 
polynomial in V(t), which characterizes the ionic current. The 
second zero is Vthr- Beyond this threshold the total ionic current 
switches its sign. So eqn. (1 1) becomes: 

IiONo(V)=-0-04V 2 -5V + bV r -\40 

(14) 

= -0.04(V-V r )(V-V THR ) 



Hence, V r = — 70 mV and the resting threshold is Vthr,o = — 
55 mV. 

We will utilize this simple nonlinear model to complete the 
picture. If U'( V ,f) > w r - i.e. the membrane is not at rest, the point 
where the total ionic current Iion{V) switches sign is shifted 
rightward toward a higher Vjhr value. For example, for very long 
durations 7"— >oo, w—*bV: 



lion,*, ( V) = - 0.04 V 2 + (b - 5) V - 1 40 
= -Q.04(V-V r )(V-V THR ) 

The subscript oo indicates that we have assumed w(t) = bV(t). 
Predictably, this does not affect the resting potential, since 
Ii<m,tx>(V r ) = Ii m fl(Vr). However, Vthr,™ = -50 mV is higher 
than the resting threshold VjHRfi- 

This reflects the lowering of excitability shortly after an AP, and 
once the post-AP membrane re-polarization takes place. This is 
known as refractoriness, which can be either absolute - i.e. no AP can 
be elicited regardless of how large the stimulation, or relative - i.e. 
larger stimulation current is required - to reach a higher threshold 
Vthr- 

Some models of the HH-type have even more complex 
Iion,oo(V) and thence Vthr behavior. This complexity is due to 
the multiple gate states, which may have very different time constants 
and hence reach their asymptotic states at different times. In 
addition, the HH models involve inactivating sodium (Na + ) 
channels. Hence, excitability may be conditional on attaining 
the firing threshold within a specific time window. Then Vthr may 
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exist only with durations « go . Hence, even over arbitrarily long 
duration, an arbitrarily low (non-zero) current may never elicit AP's, 
and may also damage the tissues and the electrodes as irreversible 
chemical reactions take place. 

So, wide stimulation pulses lasting well over some critical 
duration Tqr may not be able to elicit any AP. This is due to the 
comparable temporal scales of duration Tstim an d the time 
constant T IO „ of the closing gates associated with depolarizing ionic 
currents and of the opening gates associated with re-polarizing 
currents. 

Therefore, let us assume that the excitable-membrane's 
potential is at its resting value V r . Hence, in principle an action 
potential (AP) can be elicited by stimulation of the fixed duration 
T <Tcr. Therefore stimulation takes place over a finite time- 
horizon. 

Finite-Horizon Optimal-Control (FHOC) 

In this approach, the current waveform is the unknown system 
input signal complying with specific optimality criteria. The 
optimal pattern u*(t) for ?e[0,T] is sought as a solution of the 
following constrained minimization problem: 

u'= argminj<9(X(r))+ J f 0 (X,u)dt^ (16) 



dt 



X = F(X,u) 



Vu(t)e[L,R] 



where L and R are the constant lower and upper bounds on the 
values for each u(t) sought. 

The computational model's dynamical system is introduced in 
the optimization problem of eqn. (16) in the form of a set of 
equality constraints. The vector function F(x,w)eR" describes the 
dynamics of the array of system state-variable trajectories 
Xi(t),i = 1 . . .«, resulting from given initial state X(0) and control 
signal u. 

The example developed in the Results section uses the 
Izhikevich model - eqns. (6) and (1 1) - with n = 2. 

The minimized functional, contains the integration term fo(X,u) 
and a final-time (also known as penalty) term 0(X(7")) - pulling 
toward the desired final state X*(7"). The specific fo expression 
yields minimum electric stimulation power: 



vector-matrix transpose operator. H =fo(x,u) + A'F(X,m) is known 
as the Hamiltonian. 

The necessary conditions for optimality require that all partial 
derivatives of the Lagrangian by the system states vanish at the 
optimal solution to the problem of eqn. (16) - i.e.: 



8C 
dX(t) 



= 0 



\/te[0,T] 



(20) 



Here the 'vector-matrix' notations 8(j)/8X or 3F/3X, where 
XeR", mean respectively d(j>/8xj or dfi/dxj, Vij= 1 . . . n. 

This development is known as mathematical sensitivity analysis 
and its main purpose is to reveal the impact of a given system 
parameter (such as u(t) or its initial state X(0)) on the resulting 
dynamics. 

From eqns. (19) and (20): 



d _ dH 
Jt X ~~~SX 



(21) 



A(?> 



Be 

dX(T) 



where 



dH _df 0 5F' 
5X ~ dX + dX X 



Notice that eqn. (21) describes the adjoint dynamic system 
iterated in reverse time with a terminal condition provided by the 
derivative of the 6(X(T)) term. To solve the ODE system of eqn. 
(21), the achieved forward dynamics of eqn. (16) needs to be 
already computed. 

Similarly, all partial derivatives of the Lagrangian by 
u(t),\/te[0,T] vanish at the optimal solution to the problem of 
eqn. (16) - i.e. \fk = 0 . . .m— 1: 



dC_ 
But 



-I 



8H 



dt 



where At is the sampling time, W/t = u(kAt) and 



(22) 



MX,u) = u(t) 2 /2 



(17) 



The penalty term is a convenient way to express the desirable 
stimulation's outcome - the membrane voltage reaching some pre- 
defined threshold-level Vthr' 



6(x(7> 



K, 



penalty 



2(Vthr — V(T)) 



(18) 



Using a general constrained parametric optimal-control ap- 
proach (e.g. [32]), the objective and equality constraints in eqn. 
(16) are combined into the Lagrangian: 



du du du 

Hence, eqn. (22) yields all components of the gradient w.r.t. 
u(kAt), which enables the use of gradient-based quasi-Newton 
search routines (e.g. fmincon from the Madab optimization 
toolbox). 

Moreover, one sees from eqn. (19) that the array \(Q) is the 
sensitivity (i.e. the gradient) w.r.t. initial state X(0), i.e.: 



A(0) = 



dC 

ax(0) 



£ = 9(x(r)) + 



f 0 (X,u)dt-X(-rX-F(X,u)) 
dt 



=6(x(r)) + [A'X]f =0 + 



d , 
tf+^A'X 



(19) 



where \(t) are the Lagrange multipliers, associated to each of the 
n equality constraints in eqn. (16) and (•)' stands for the 



A boundary-value problem (BVP), with known initial conditions 
for X(0) and terminal conditions for 1(T), is solved numerically. 
However, it should also be noted that such solutions may also 
converge to shallow local minima. For example, the Newton search 
is guaranteed to produce the 'true' solution when the problem at 
hand involves a quadratic cost. Here the objective function not 
only may be non-quadratic, but also may be non-convex in some 
manifolds of its high-dimensional parametric space. 
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Figure 3. Propagating AP's and spatial profile of the membrane voltage V(t,z) & intracellular potential <\>J Isi iai (at the end of 
stimulation, please also see Fig. 2); z is the 1 D axonal spatial coordinate. The peaks of V at the Ranvier nodes are due to the direct exposure 
to the extracellular medium, which is unlike that of the myelinated sections in the double-cable MRG'02 model. 
doi:1 0.1 371 /journal.pone.0090480.g003 



Above we described the continuous-time FHOC. The CVODES 
toolbox readily provides adjoint sensitivity analysis (ASA) capabil- 
ities. FHOC is one of the common applications of the latter. 
Analogously, a discrete-time version may be formulated and solved 
(see the Results section, where a specific example is developed). 

Solving the Problem Analytically: The PLA in E5 

Through calculus of variations, here we establish a general form 
for the energy-optimal current waveform u*(i). This approach 
applies the Principle of Least Action to ES. 



Let us assume that T«Zjon, where Tion is the time-constant 
that determines the behavior of the slow gate states of the modeled 
ionic-channels. Hence, the fast gate states may be approximated by 
their asymptotic values x co (V)= lira t^ m x(t\V), while the slow 
gate states - by their resting values Xo=X aj (V r ). 

Then an AP can readily be evoked by stimulation from the 
resting state, and the threshold potential Vjhr to reach at time T 
is finite and assumed (without loss of generality) to be known. The 
energy-efficiency of driving the excitable-tissue membrane poten- 
tial V(t) from its resting value V r to Vthr through a stimulation of 
fixed duration T satisfies: 
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Table 6. MRG'02 double-cable model-axon geometric parameters, in \.im. 





Notation Parameter description 


Value 


Shared parameters: 


D Fiber Diameter 


16.0 


AZ Node-node separation 


1500 


N my Number of myelin lamellae 


150 


Nodal compartments: 


L„ Node length 


1.0 


d l7 Node diameter 


5.5 


MYSA (myelin attachment paranode) 


Lm length 


3.0 


d M diameter 


5.5 


S M periaxonal width (Membrane-to-Myelin gap) 


0.004 


FLUT compartments (main section of paranode) 


L F length 


60.0 


dp diameter 


12.7 


Sp periaxonal width 


0.004 


STIN compartments (internodal section, 3+3 total in 1 internode, see Fig. 2) 


L s length 


228.8 (*1 


ds diameter 


12.7 


S s periaxonal width 


0.004 



doi:1 0.1 371 /journal.pone.0090480.t006 



u*(t)= argminP(«) P{u)=\/2 



[u{t)] 2 dt (23) From eqns. (24) and (26), and since u\i)=C m V* +h(V*). 



Since from eqn. (6), u(t) = C m V + I-z(V): 



P(u) = S(V\T,u)= 1/2 



[C m V(t)+h(V)] dt (24) 



As done in the calculus of variations let us perturb the energy- 
optimal time-course K*(Z)bythe infinitesimal perturbation Ol(t), 
where r\(i) is an arbitrary function of time and e is an infinitesimal 
scalar. 



P(e) = S(V*)+e 



-<?F(V*,tj) 



u*(t)(C m i,+tiI^V*))dt 



(27) 



The necessary condition for S( V) to have a minimum at e = 0 for 
any ;;(/) is: 



G e = P'(e)| e=0 = u"(,t)(C m fi+r,I^(V*))dt=0 (28) 
Jo 



V(t)=V*(t)+er,(t) 

dIy(V*) 

oV 



(25) 



From eqn. (25), \/te[0,T] the integrand in eqn. (24) becomes: 

(C m V+h( V)f = (C m V* + h( V*)) 2 

+ 2e(C,„V* +h(V t ))(C m i 1 + , 1 I^ (V)) (26) 
W(C m r, + r,Ii ( V*)) 2 



To deal with the u*(t)rj term of eqn. (28), it is integrated by parts : 



G e =c m [u*(t)n(t)]l 



[C m u* -u*l! L (V*)]n(t)dt = Q (29) 



Since the perturbation r\{t) respects the boundary-value 
problem (BVP) with known initial and terminal conditions for 
V*(t) - i.e. r](0) = ri(T) = 0, then the first RHS term above 
vanishes. Hence, the only way that eqn. (29) will hold for any rj(t) is 
that we have the Euler-Lagrange-type equation: 
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Figure 4. LAP energy-optimal I ! w) and u*(t) for the LM: for T S tim respectively 10 [is and 5 ms; the time constant z=C/g was varied as 
indicated in the legend; membrane capacity was constant - C=1 fiF/cm 1 , while membrane (leak) conductance g was respectively 
0.2, 1 and 5 mS/cm 2 ; The 3 solutions shown correspond to the nominal i = 1 7ns (cyan trace) or 5-fold shorter (thin red dash-dot), or 
5-fold longer (thick dashed black) t respectively; (thin dashed black) rectangular pulse with amplitude k*= {V THR — V,)/T S tim- 
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c m u* = i%(y*)u* 



(30) 



H=J/2+\{u-h{Vj)IC m 



(31) 



Equation (30) can also be attained directly using the continuous 
version of the standard OC formalism [32] (please see also the just 
presented FHOC subsection above). 

Here the Hamiltonian is. 



The necessary conditions for optimality require that. 

dH/du = 0 



(32) 
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Figure 5. LAP optimal waveforms V*(t) and u*(t) for the OD IM: 
The 3 solutions shown correspond to the nominal IM opposing 
current (cyan trace), twice higher (thin red dash-dot), or twice 
lower (thick dashed black) /v respectively. The Iwn,q{V) 

approximation of the ionic current is used for a case of very short 
duration (Tstim = 10 /is) and the Iion,w( V) approximation is used for a 
case of long duration (Tstim = 5 ms). It is important to notice that - as 
with the LM model above, u*(t)^k* + Iion(V), where 
k* =(V THR -V r )/T (see the Box) Box: Resting-state I IO n,o( v ) and 
asymptotic-state Iion,x{V) ionic currents for the 0D IM; Markers are 
inserted at the resting and threshold membrane-voltage points, 
respectively Vrest = ~70, Vthr,o = ~ 55 and Vthr,x, = ~ 50 mV. 
doi:1 0.1 371 /journal.pone.0090480.g005 



C 2 m V*=I^(V*)[u*-C m V*}. 
And finally, from eqn. (6). 



^=/ S (nx^ (34) 

Equation (34) is a rather simple system of ordinary differential 
equations (ODE) that can readily be solved for a given current 
model Iz(V*) to compute the energy-optimal membrane voltage 
profile V*(t). The energy-efficient current waveform u*(t) is then 
computed from eqn. (6). 

In the Results section below we illustrate the use of eqn. (34) 
with several frequently encountered current models. 

Results 

Here, we first derive some key analytical results using the 
simplest and clearest models. We then identify generally applicable 
optimality principles. Finally, we demonstrate how these principles 
apply also to more complex and realistic models and their 
simulations. 

Part I - Specific Point-model Results, Applying the LAP 

For the zero-dimensional (single-compartment, space clamp) 
models introduced in the Methods, here we describe the LAP- 
optimal waveforms V*(t) and u*(t), stemming from the general 
(model-independent) LAP result of eqn. (34). 

These simple cases readily illustrate some rather key optimality 
principles resulting from a LAP perspective. We will discuss these 
optimality principles as we go, and will summarize them at the end 
of this subsection. 

Linear sub-threshold model. Replacing Iy.(V*) in eqn. (34) 
with Iion(V) from eqn. (8): 

x 2 V=V (35) 

T = C m /g m = R,„C m is the membrane's time constant and for 
expediency V =V* and V r = Q. 

The general solution of eqn. (35) is: 



\= -8H/8V 



(33) 



From eqns. (32) and (31) \/C m = — u. Then from eqn. (33). 



C m ii = dH/dV= -\l C m U V) = ul'zi V) 



which is the same as eqn. (30). 

From eqns. (6) and (30) we have that. 



U*(t) = C m V*+I^(V*)V* =I'^V*)u*/C m 



and thence: 



V(f) = C l e*+C 2 er 



(36) 



Given the boundary conditions V(0) = 0 and V(T)= Vthr 

sinh(r/T) 



V*(t)= V T 



' sinh(r/t) 



(37) 



A result similar to eqn. (37) is obtained by [33], using a slightly 
different (less direct or general) optimal-control approach. 

From eqn. (37) one can see that V*{t)/ Vthr = 
sinh (t/x)/ sinh (T/t)^t/T - i.e. it has a linear rise, especially with 
r«T. Here T= 100 fis and t=1 ms (computed using typical 
values from the literature for g m =\ mS/cm 2 and C„, = 1 
fiF /cm 2 ). 
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Figure 6. LAP optimal waveforms V*(i) and »' (/) for the OD HHM: The Iion,o( 1 > approximation of the ionic current is used for a case 
of very short duration (T S tim = 10 fis) and the IroN,w{V) approximation is used for a case of long duration (T S tim = 5 '«■') (see the Box). 
As with the IM, bvp4c was used to numerically solve the BVP of eqn. (34). The figure follows a quite similar format to Fig. 5. /^(K) can also be 
assumed higher or lower. All the maximal ionic conductances in the HHM (see also Table 3) are temperature-dependent and are linearly proportional 
to the coefficient kj. The 3 solutions shown correspond to the ionic current at T°C\¥I°C\ (cyan trace), twice higher (thin red dash-dot), or twice 
lower (thick dashed black) h respectively. From eqn. (42) we can see that kr = 1.6047 (half the nominal) at 28.7° C, and /(r = 6.41 88 (twice the 
nominal) for at 45. 3° C. Box: Resting-state 1ion,q(V) and asymptotic-state Iion,co(V) ionic currents for the OD HHM; Markers are inserted at the 
resting and threshold membrane-voltage points, respectively Vrest = —llmV, Vthr.o = ~ 64.55 mV and Vthr.oj = ~ 52.35 mV. 
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V(i) from V r to Vthr at a constant rafe, is the time-constant 
waveform u(t) = k* = (Vthr — V r )/T. For this example, k* = 
\Q»Iion(V), which explains why u*(t) is that close to a 
rectangular waveform. 

As a matter of fact, for very short stimulation times, the k* tend 
to be high, while Iion(V) tends to be linear. Hence, the 'classic' 
rectangular (or square, SQR) waveform tends to also be close to 
energy-optimal. 

Such facts are rather important as they lead us below (as 
evidence is accumulated) to a general form not only of V*(t), but 
also of u*(t). 

Comparative properties the V(t) growth profiles. The 

GE waveform may be an SQR waveform in disguise. I.e. some 
linear growth of the membrane voltage may still fit the one 
obtained upon ES with a GE. The motivation for this is in eqn. 
(36), where the first term vanishes with T»x. 

Finally, the total electric charge conveyed by the ES source may 
have to be considered. For example, in the LM of eqn. (8) the total 
charge consists of a capacitive charge to raise the membrane 
voltage by a given amount (to Vthr), and resistive charge 
jTsrm V (t)jRdt. A similar situation occurs in the MRG' 02 model 
due to the opposing axial currents. 

So let us solve the following auxiliary problem: 
Find a linear fit V(t) = max[a(7 — b),0] to the growing exponent 
K(?) = (e' /l -l)/(e 7 '- s ™ /T -l), so that the ES source conveys the 
same resistive charge in the time interval te[0,TsTm]- Le. we want 
that: 



-80 -60 -40 -20 0 20 40 60 
V [mV] 



V(t)dt-- 



V(t)dt-- 



Tst 



e T STIM/ T — l 



Figure 7. The LAP vs or with numerical optimisation for the OD 
IM, with I siim = 2 ms: see also Fig. 5 which shows that an initial 
guess u*(t), based on the //near-growth rate 

k* =(Vthr— K)/Tstim is still valid with T S tim = 2 ms and 
Vthr = -50 mV. panel A: discrete-time IM and FHOC panel B: 
continuous-time IM and FHOC, using CVODES adjoint sensitivity analysis 
capabilities upper plots: (dashed black) a rectangular pulse with 
amplitude k*\ (thick cyan) the LAP — k* -\- Iion.oi 

(V); (thick black) 

the best FHOC u(t) lower plots: (dashed black) //near-growth evolution 
of the membrane potential from V r at t = 0 to Vthr at t=TsnM', 
(dotted gray) the desired threshold value Vthr = ~ 50 mV; (thick cyan) 
the resulting LAP V*(t); (thick black) the resulting FHOC V(t). 
doi:1 0.1 371 /journal.pone.0090480.g007 

Figure 4 presents the LAP energy-optimal stimulation profiles 
V*(t) and u*(t) for a short and a long stimulus duration Tstim 
and three membrane time constant x values. 

Before we go on, it is useful to investigate the conditions for a 
growing exponent (GE) waveform to outperform the SQR 
waveform. 

First, u* GE (t) has a very rapid rise. Hence, its optimal duration 
T GE will be short. Second, it is noteworthy that in [33] x — 30.4 
micro-seconds! Hence, injected current rapidly leaks out. However 
even with the above extreme x value, at its optimal duration T* S q R 
the SQR wave does just 22% worse, which means that the SQR is 
among the best candidates for its robustly good performance. 

Second, in multiple cases, the energy-optimal LAP waveform 
u*(i) looks a lot like a 'classical' rectangular waveform. From eqn. 
(8), we may also see that, with V r = 0, Vthr = 1 , the max. value of 
Iion( V) is equal to 1 and is attained as the membrane potential 
reaches the threshold V(T)=Vthr- If we then replace 
Iion( V)x0 in eqn. (6), we see that a waveform u(t) - that brings 



Here, for simplicity (and without any loss of generality) we have 
assumed V r = 0 and Vthr = 1 ■ 

For example with x= Tstim /4, we obtain 6«0.54 x Tstim, i-e. 
the linear-growth equivalent has more than twice shorter duration 
- e.g. with Tstim = 5, A«2.7. 

The latter result promotes intuition: with large opposing 
currents optimal ES cannot afford to last long. The transition of 
the membrane voltage from its rest to a threshold value is best 
performed rapidly. Hence, the shape of the V(t) growth profile 
depend on the Tstim I x ratio. As seen, for Tstim «x, the optimal 
u(t) is close to rectangular, while with Tstim »x, the GE is in 
effect equivalent to doing nothing for at least half of the duration, 
and then to a SQR waveform of at least doubled amplitude. 

With quite similar reasoning, one can demonstrate that a lst- 
order membrane voltage growth profile V(t) = {\ — e'' 1 )/ 
(1— e Ts ™/ T ) in the time interval te[0,TsTiM] is suboptimal and 
equivalent to linear growth, which has about twice longer 
duration. 

Izhikevich model. Replacing in eqn. (34) with the 

Iion(V) approximations from eqn. (14) or (15), see Box in Fig. 5: 

C 2 m V = QM 2 (V- V r \V- V THR )[2V-(V r + Vthr)} (38) 

As in the preceding model V=V*. Note that the dynamics of 
eqn. (38) has all FP's of Iion(V), as well as a third FP at 
V = 0.5(V r + Vthr), contributed by the derivative term I' ION (V). 

Equation (38) can be solved analytically. However, it provides 
the solution in an implicit form and involves an incomplete elliptic 
integral of the first kind. Hence, we used the Matlab bvp4c BVP 
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Figure 8. The MRG'02 model: Toward upper bounds on V thr (T S tim)' the figure presents a family of ionic current Iion(V,Z) 
approximations at the target site (Z = 0), for a set of durations T S tim- For each of the durations it is assumed that the membrane voltage 
trajectory V(t) evolves according to a linear ramp from rest V r to threshold Vthr (the unknown). For each V value on the horizontal (independent- 
variable) axis of the figure, a V(i)=kt ramp was assumed and the corresponding ionic current Iion(V) was computed, based on approximate gate 
states (see the Box). Note: for the sake of better visibility, a x 10 gain is applied to the approx. 1ion(V) for the case of Tstim = 5 ms. Box: For a 
chosen Tstim = 5 ms and as V(t) is linearly ramped up, for each gate state the plots show the ratio x( V, T)jx m ( V), where x is given by eqn. (46) to its 
asymptotic value - both functions of V. Legend for gate states: opening in and closing h gates for the fast Na + ion-channel subtype; p persistent Na + 
channel gates; s slow K + gates. 
doi:1 0.1 371 /journal.pone.0090480.g008 



solver with boundary conditions V(0)=V r and V* (Tstim) = 
Vthr- 

Figure 5 illustrates the energy-optimal LAP solution u*(t) and 
the corresponding membrane voltage profile V*(t). The IioNfli V) 
approximation of the ionic current is used for a case of very short 
duration (Tstim = 10 fis) and the Iion,co(V) approximation is 
used for a case of long duration (Tstim = 5 ms). 

It is important to notice that - as with the LM model above, 
u*(t)xk*+I ION (V), where k* =(V TH r - V r )/T (see the Box in 
Fig. 5). 

According to eqns. (14) and (15) the opposing current in the IM 
can be presented in the general form: 



h(V) = gainxI I0N (V) 



(39) 



u dtr 



V L dt 



(40) 



By the Cauchy-Schwartz inequality in the space of continuous real 
functions, it is straightforward to show that the voltage trajectory 
V*(i) that minimizes eqn. (40) is such that V*(t) = k* , where k* is 
determined from the boundary conditions satisfied by V*(t). 
Hence: 



k*- 



Vthr — V r 
Tstim 



(41) 



where the nominal gain=\, and Iion(V) = 0.04(V — V r ) 
(V-Vthr)- 

To see how the optimal ES is affected by the level of opposing 
current, it is more than tempting to experiment with different gain 
values. 

Hence, 3 gain cases are plotted in Fig. 5 - for the nominal gain 
(cyan traces) and two additional cases: the opposing current 
Iz(V*) is either doubled (gain = 2, red traces) or decreased two- 
fold (gain = 1/2, black traces). As could be intuitively expected 
from the general equation (24), when Iion(V)^>0 (very low ionic 
currents): 



Just as in the preceding model, it is also V* (t) j V thr^ 1 1 T stim 
with the shorter durations - which justifies the use of the resting 
approximation IioNfi( V). 

HHM. Here the h(V*) of eqn. (34) is replaced with the 
resting-state - IioNfi(V), or asymptotic-state - IioN,m(V) ionic 
current approximations (see the Box in Fig. 6). 

Toward hoNfi(V) the gate-state variables are factored out as 
follows: The fast state wi(?)a;m 0D (F'), while the slower variables 
h(t)Kh r = hrx,(Vr), and n(t)Kn r = n ca (V r ) are approximately at 
rest, assuming very short durations. Conversely, and assuming 
very long durations, toward hoN,m{V) all gate variables are 
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Figure 9. The actually computed I ;7/; , as a function of I sum : Notice how the computed r,, ;/; value is rather similar (almost 
matched) between the linear and exponential cases, for I sum respectively 2 and 5 ms; and between the I -order and linear cases, 
for Tstim respectively 0.2 and 0.5 ms. see also Fig. 1 0. 

doi:1 0.1 371 /journal.pone.0090480.g009 



approximately at their asymptotic value, corresponding to a given 
membrane voltage V(t) (see Methods). 

As with the IM, we used bvp4c to numerically solve the BVP of 
eqn. (34) with boundary conditions V*(Q)= V r and V*(Tstim) = 
Vthr- 

Figure 6 follows a very similar format to Fig. 5. 

Similarly to eqn. (39) above, IziY) can be assumed higher 
or lower. All the maximal ionic conductances in the HHM (see 
also Table 3) are temperature-dependent and are linearly 
proportional to the coefficient kr' 

kT = Q\ T 0 - T ° )n0 (42) 



where Qio =2.3 and To = 23°C. Hence with T = 37°C, according 
to eqn. (42) k T = 3.2094. Let this be our standard case (gain = 1). 

As we did with the IM, 3 gain cases are plotted in Fig. 6 for 
Iz(V*) = gain x I [ON (V). For the two additional cases the 
opposing current is either doubled (gain = 2, red traces) 

or halved (gain= 1/2, black traces). 

Once again - as with the LM and IM models above, 
u*(f)xk* +Iion(V) (see the Box in Fig. 6). 

Numerical model simulation and optimal control. The 
IM was also evoked in the FHOC Methods section. It is therefore 
interesting to contrast the results of the LAP and FHOC 
approaches in identifying energy-optimal ES waveforms for the 
same ionic current model. For such comparison, the IM has the 



Table 7. Minimal Vnni(Tsmf)[fnV\ values for the MRG'02 model, obtained for each V(t) trajectory class. 





Tstim 


Linear 


Tst-order 


Exponent. 


0.020 


-25.649 


-37.602 


-4.963 


0.050 


-41.838 


-50.515 


-24.311 


0.100 


-50.852 


-57.366 


-37.032 


0.200 


-57.061 


-61.506 


-47.137 


0.400 


-60.588 


-63.558 


-54.124 


0.500 


-61.247 


-63.889 


-55.731 


1.000 


-62.378 


-63.960 


-59.255 


2.000 


-61.950 


-62.578 


-60.977 


5.000 


-59.273 


-59.094 


-61.249 
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Figure 10. The energy P and charge-transfer Q values as a function of Tstim ■ The linear-ramp voltage profile yields the best P 
performance for most of the durations. As in Fig. 8 notice that the P and Q values are quite similar for the linear and exponential cases, for 
Tstim respectively 2 and 5 ms; and also for the l"-order and linear cases, for Tstim respectively 0.2 and 0.5 ms. Toward the P values electrode 
impedance of 1 M 0 is assumed. Contrasted: SQR stands for the square (or rectangular) stimulation waveform. 
doi:1 0.1 371 /journal.pone.0090480.g010 



clear advantage of hiding no implementation specifics inside a 
black box. 

The FHOC formalism (see Methods) is computationally 
efficient, but it is also subject to the similar limitations as most of 
the ad-hoc search approaches. Iterative numerical optimization 
requires an initial guess for the solution, and trying different 
starting arrays 2/ 0 ' may alleviate a bit the propensity to converge to 
shallow local energy-minima. 

Here it is also important to realize that in eqn. (16) the two 
terms to minimize in the F(u) functional (a function of functions), 
namely the energy cost (17) and the penalty ( 1 8) may conflict each 
other. When the penalty gain K pena i ty in (18) is too low, the search 



will identify a lower-energy solution u, which however does not 
bring the membrane potential Vt up to the desired threshold value 
- i.e. Vm « Vthr- Conversely, a too high penalty gain K pma i ty will 
identify a very high-energy solution it, which is not only costiy, but 
the membrane potential may also overshoot the threshold, since 
the 'getting there' is underestimated for the sake of the very last 
simulation steps. 

As seen from Fig. 7 Panel B (which uses the Iion,o>(V) 
approximation of the ionic current for the relatively long duration 
Tstim = 2 ms), the linear growth profile is a reasonable estimate 
for the optimal membrane voltage profile V*(t). Hence: 
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Table 8. Minimal Q(T ST!M )\pico- 


Coulomb] values for the MRG'02 model, obtained for 


each V(t) trajectory class. 




TsTIM 


SQR 


Linear 


1 st-order 


Exponent. 


0.0200 


3.1180 


0.1279 


0.1671 


0.1467 


0.0500 


1 .2472 


0.1630 


0.1946 


0.1642 


0.1000 


0.6236 


0.1959 


0.2212 


0.1847 


0.2000 


0.3832 


0.2369 


0.2583 


0.2121 


0.4000 


0.3426 


0.3045 


0.3191 


0.2545 


0.5000 


0.2605 


0.3440 


0.3492 


0.2937 


1 .0000 


0.2143 


0.5093 


0.4736 


0.3910 


2.0000 


0.1808 


0.8640 


0.6361 


0.5855 


5.0000 


0.1411 


2.1216 


1.5673 


1.2018 


doi:l 0.1 371 /journal.pone.0090480.t008 



u*(t)xk*+I h 



(43) 



where k* is given by eqn. (41). When w' 0 ' is close to the LAP 
estimate u*(t) of eqn. (43), the FHOC iteration also consistently 
ends close to there (see Fig. 7, panel B). The cyan traces on Fig. 7 
are the u*{t) and the resulting V*(t). With the LAP estimate, the 
FHOC approach resulted in a final membrane potential 
reasonably close to the desired threshold value - i.e. 
V(Tstim) = — 50.106 w Vthr = — 50, even if the IM was simu- 
lated with the discretized LAP waveform u*(t) (Af = 10 fis). 

The black traces illustrate the FHOC solution, computed for 
two different i/ 0 ' choices. For Panel A, was chosen to be all 
zeros. When all time-step entries M* 0 ' were chosen to be equal to 
the upper bound U = 30 (data not shown), due to the (discontin- 
uous) AP event occurring mid-way the temporal horizon, the 
Matlab's fmincon solver remains stuck to the initially provided 
values. 

Except for the case in Panel B, the K pena i ty meta-parameter had 
to be kept high [Kp em i ty = 70) in order to respect the terminal 
constraint of V(Tstim)~ Vjhr- 

The total energy costs (all expressed as 2-norms of the obtained 
best u(t)) are respectively 161, 153.2 and 423.4 (for the discrete- 
time version) 186.7, 159.1 and 334.2 (for the continuous-time 
version). 

Comparing these to P(u*)= 153.2 (discrete-time) and =157.4 
(continuous-time), the LAP-based solution is comparable to or 
superior than the FHOC solutions. The numerical FHOC solution 
on Fig. 7, panel A has converged to a local extremum. Note that a 
post-hoc correction (simple DC offset) is applied to the LAP-based 
estimate, which adjusts for the overshoot of Vjhr when simulating 
the full (two-dimensional) IM. The overshoot is due to using the 
one-dimensional approximation, eqn. (15). 

The results obtained here nicely illustrate multiple aspects of 
identifying energy-efficient waveforms through numerical model 
simulation and optimization. Clearly, pairing theoretical insights 
with numerical tools carries the best success potential. 



good idea of both Iz( V) and Vjhr to successfully solve for V*{t), 
and thence for u*(t), in a given particular situation. 

We identify also the following key and practice-oriented 
optimality principles resulting from the LAP perspective. 

1 . The optimal sub-threshold membrane potential growth profile 
with relatively short durations Tstim and low membrane 
conductivity: 

First, in all simple models we used up to here, the solution V*(t) 
of the ODE system, defined by eqn. (34), is quite close to a linear 
growth from V*(0) = V r to V*{Tstim) = Vthr- Second, with the 
total current « 0 (e.g. low leak), then from eqn. (6), it 

follows that u(t) will be exactly proportional to the rate of 
change of the membrane's potential V(t). If V*(f) « const, then 
u*(t) is close to a SQR waveform. 

2. The energy-efficient waveform depends directly on the 
temporal shape of currents at the AP initiation site. 

3. The targeted Vjhr membrane voltage threshold depends on 
stimulation duration, with a tendency to increase with Tstim- 

4. The exponential growth membrane voltage profiles V(t) are 
equivalent to linear growths of shorter duration. 

Part II - Multiple-compartment Model Results 

Here we first extend the general (model-independent) LAP 
result of eqn. (34) to spatial-structure models (non-zero-dimen- 
sional, multi-compartment), which involve membrane-voltage 
distribution and propagation along cable structures. 

LAP result generalization to multi-compartment 
models. There is a combinatorial explosion in both the number 
of parameters and the number of ways that multi-compartment 
models can be put together and used. Hence, there is much more 
than one way of generalizing the LAP result of eqn. (34). 

Here we briefly present a variant, which appears to be one of 
the most straightforward generalizations. 

With a multi-compartment model, eqn. (7) can be rewritten as: 



Part I Results Summary 

A number of more general observations on u*(t) can be made 
looking at the results this far. 

Probably, the most significant result is that the use of LAP 
reduces the problem to the BVP, defined by eqn. (34), with 
V*(0) = V, and V*(Tstim) = Vthr- We still need to have a very 



jV(t,Z) = u(t,Z)-h(V,Z) (44) 

Without loss of generality, we used the variable Z to represent 
any 'spatial' model dimension. It could even stand for the 
compartment index in a discretized implementation. 
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Figure 11. Optimal waveforms u*(t), Tstim = 20, 200 us: The 
figure also provides the corresponding optimal SOR-\\ke 
linear-growth-related current C„, x k (dashed black), as well 
as the components of /v - respectively the 1 ]0 n (blue traces) 
and /,,,,,; (red traces) current trajectories. 
doi:1 0.1 371 /journal. pone.0090480.g01 1 

Now, eqn. (7) is a partial DE, depending both on the temporal 
and the spatial model dimensions. 

Assuming that we are free to manipulate u(t,Z) in every 
compartment as we wish, the derivation sequence from eqn. (23) to 
eqn. (30) (see the LAP subsection in the Methods) still applies 
yielding a family of equations 'parameterized' by the location 
coordinate Z. 

Hence, we may obtain the generalization of eqn. (34) as: 

C 2 m ^V*(t,Z) = h(V\Z)x^h(V\Z) (45) 

Like the extended eqn. (44), eqn. (45) is a partial DE, depending 
on both temporal and spatial boundary conditions. In particular, 
Vthr becomes a function of Z. It is no longer a single variable, 
but a whole spatial profile, subject to conditions such as the safety 
factor for propagation introduced in the cardiac literature [34] . 

The MRG'02 model: Toward upper bounds on 
Vthr(Tstim)- Multi-compartment models add complexity 




t[ms] 

Figure 12. Optimal waveforms »'(/): see also Fig. 11. Notes: 

Since here Vt,V*'(t) = k t , where k* is given by eqn. (41), from eqn. (6) 
u*(f) = C m k* +h(V). The figure is optimized to present clearly both 
u*(t) and k" (*1) The dashed trace at the bottom plots log 10 C m k* as a 
function of Tstim (*2) Toward equally good plot visibility, for all 
durations Tstim < \ms, the waveforms u*(t) are rubber-banded to take 
the same graph width as the 1 ms-waveform. This is illustrated by the 
scale bars for the shortest duration Tstim = 20 /i.?. (*3) The vertical scale 
is the same for all plots, except for the logarithmic offset, as defined by 
pt. (*1) above. 

doi:10.1 371/journal.pone.0090480.g01 2 

unseen with the single-compartment models. Wongsarnpigoon & 
Grill [8] used the peripheral-axon MRG'02 model [25] in a 
genetic-programming search for energy-efficient stimulation 
waveforms. The approach was somewhat similar to the FHOC 
described above. After thousands of iterations simulating the 
MRG'02 model, the identified waveforms were reminiscent of 
noisy truncated and vertically offset Gaussian's (Fig. 2 in [8]). In 
the light of analysis this far one might think that this reflects the 
shape of V) for V ranging from the resting value (—80 mV) to 
some threshold Vthr- 

In this work stimulation is assumed to be intracellular and at just 
one spatial location (Z = 0, the center RN, see Methods) along the 
cable structure. 

To suggest a version of optimal waveforms u*{t) for the 
MRG'02 model, we first estimate the membrane voltage threshold 
for each duration. One analytic way toward such estimates is 
readily provided by the MRG'02 model. Recall also that with 
simpler models Vthr showed a tendency to increase with Tstim- 

Figure 8 presents a family of ionic current Iion(V,Z) 
approximations at the target site (Z = 0), for a set of durations 
Tstim ■ For each of the durations we assume that the membrane 
voltage trajectory V(t) evolves according to a linear ramp from 
rest V r to threshold Vthr- As the latter is unknown, we produced 
one such ramp for each V value on the horizontal (independent- 
variable) axis of the figure, and then computed the corresponding 
ionic current Iion(V) as described next. 

Toward gross estimates of Vthr, we first solve approximately 
eqn. (10) for each gate-state: 

x( T) = x 0 + (x m ( V) - *o)( 1 - e - ™ V) ) (46) 
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Table 9. Minimal P(T STm )\femto- 


Watt] values for the MRG'02 model, obtained for 


each V(t) trajectory class. 




TsTIM 


SQR 


Linear 


1 st-order 


Exponent. 


0.0200 


1.9444 


0.9620 


1 .4387 


2.2993 


0.0500 


0.7778 


0.6391 


0.7765 


1.1611 


0.1000 


0.3889 


0.4596 


0.5158 


0.7325 


0.2000 


0.2937 


0.3307 


0.3692 


0.4766 


0.4000 


0.2934 


0.2693 


0.3003 


0.3352 


0.5000 


0.3392 


0.2675 


0.2913 


0.3463 


1 .0000 


0.4593 


0.2934 


0.2929 


0.2954 


2.0000 


0.6535 


0.4265 


0.3321 


0.3204 


5.0000 


0.9949 


1.0339 


0.9486 


0.5263 


doi:1 0.1 371 /journal.pone.0090480.t009 



where Xrj is the gate-state value at rest and V=(V r + V)/2 is the 
average excursion from the resting membrane voltage. 

Figure 8 shows the obtained approximate ionic currents 
Il01f(V(t)) as a function of just V for three very different 
durations - T S tim - 0.02, 0.5 and 5 ms. For T S tim = 5 ms, the 
Box in the same figure illustrates the estimated proportions-to-rest 
x(Tstim\ V)/x co (V r ) for each of the 4 gate-state variables, at the 
end of stimulation. 

Why does such an analysis provide upper bounds on 
Vthr(Tstim)? 

First, from the Box of Fig. 8 we can see that indeed the 
dynamics of the fast Na + ion channel subtype evolves before that 
of the other ion channels. Particularly, we see that the estimate for 
inactivating h gates suggests they are completely closed for 
Tstim = 5 ms and once V reaches around — 40 m V. 



T ST,M= 100 ^ 




0.2 



0.4 



0.6 
t [ms] 



0.8 



Figure 1 3. Propagating AP due to an optimal SQR (rectangular) 
waveform, Tstim = 100 (is: For the shortest durations, the plain 
rectangular waveform outperforms by P the ones associated to 
the linear-ramp voltage profile. One can see clearly that the steep 
rise of the SQR waveform yields an early superlinear ramping of the 
membrane voltage. However, the rectangular waveform requires a lot 
more charge Q to be transferred (see Fig. 10). 
doi:1 0.1 371 /journal.pone.0090480.g01 3 



On the other hand from the main Fig. 8, one can see that this 
analysis gives the intervals Ve[V re st, Vvb] in which the approxi- 
mate ionic currents Ijon(V)<0 (i.e. remain depolarizing). 

Clearly if Vthr(Tstim) is not reasonably within ]V rest , Vub], no 
miracle would yield an AP at the target location, since Iion 
becomes repolarizing outside of these bounds. 

Interestingly, the analysis also predicts lowering of Vthr with 
longer durations. This result is exactly the opposite of what was 
observed with the simpler models of the HH-type, where Iion was 
repolarizing for Ve ] V rest , Vthr] ■ 

The numerical experiments we conducted were fully consistent 
with the above predictions, and some upper bounds were also 
quite tight. 

The MRG'02 model: numerical experiments. We con- 
ducted four series of numerical experiments in search of the 
optimal waveforms u*(f) for the MRG'02 model. Each series was 
computed for the same set of 9 durations Tstim = 20, 50, 100, 
200, 400 and 500 fis; 1,2 and 5 ms (for the sake of better visibility, 
only the most representative subsets are illustrated in full detail). 

The four series differed by the chosen voltage-clamp temporal 
growth profile V(t,0) at the targeted RN location and A baseline 
series involved finding the threshold rectangular stimulation 
amplitude. In all series, the constraint was to observe a 
propagating AP at the latest within 1 ms after the end of 
stimulation. 

With AV = Vjhr(Tstim)— V„ where the minimum Vthr was 
found (with 0.001 mV tolerance) using the same type of golden- 
section search algorithm as per the optimal SQR amplitude. 

And the three LAP-driven series were: 
linear growth 



V(t)=V r + AVt/ Tstim 



(47) 



exponential growth 

V(t) = V r +A IV /T - \)/(e T sTiMh _ i) (48) 

1-st order growth 

V(t)= V r + AV(\-e-' lx )/{\-e- T STiMh) (49) 
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The corresponding u(t,0) ES waveforms were computed from 
eqn. (44) with Z = 0. 

The MRG'02 model: numerical results. Figure 9 and 
table 7 illustrate the obtained Vthr as a function of Tstim- 

The computed optimal values of Vthr ar e often similar for two 
adjacent durations either between the linear and 1-st order, or 
between the linear and exponential growth (EG). 1-st order is 
usually similar to its right-hand linear neighbor (for the next longer 
duration). Conversely, EG is similar to its left-hand linear neighbor 
(for the previous shorter duration). 

This is consistent with and best interpreted in the light of our 
growth-profiles comparison (see the dedicated subsection on page 
13). There we saw that indeed an EG V(t) trajectory is 
approximately equivalent to linear growth of about twice shorter 
duration. As for 1-st order growth, clamping the voltage to its 
plateau will tend to be similar to a linear growth of about twice 
longer duration. Recall also that 1-st order is the 'reverse-time' 
analog of EG. 

Figure 10 and tables 7, 8 illustrate the obtained optimal- 
waveforms' energy P and charge-transfer Q values as a function of 
Tstim- 

The linear-growth strategy is the one that tends to perform best 
across the board, except for the 2 longest durations, and as 
predicted by the comparative (linear vs exponential growth) 
analysis, based on the 0D LM. 

Figure 3 illustrates the propagating AP's, corresponding to the 
two representative linear and exponential voltage-clamp temporal 
growth profiles at the stimulation site V(t,0). The figure also shows 
the spatial profiles of the membrane voltage and intracellular 
potential at the end of stimulation for the two growth cases. 

Consistently with the analysis in the subsection on the 
comparative properties of the V(t) growth profiles, we found out 
that the spatial distributions of membrane voltage and intracellular 
potentials at the end of stimulation were reasonably similar - e.g. 
between the optimal linear growth voltage-clamp for Tstim = 2 
ms, Fig. 3 (Panels A, C) and the optimal exponential growth with 
Tstim = 5 ms, Fig. 3 (Panels B, D). 

Note that we expect from an approximately globally optimal 
stimulation waveform u*(t) to yield a specific distribution of 
membrane voltages V{Tstim,Z) at the end of the stimulation. 
We call this distribution tentatively the invariant spatial profile of 
the membrane voltage. Importantly, such a profile will differ for 
any different duration Tstim even when the corresponding 
waveform u*(t) is globally optimal. This is due for example to 
the small spatial constant X, which controls the spatial diffusion 
with time. 

However, if the spatial profile is about the same for different 
durations Tstim an d the corresponding different waveforms u*(t) 
(see Panels B and D in Fig. 3), then both waveforms may be 
optimal. Recall that linear fits to both the optimal 1-st order 
growth and the optimal exponential growth with durations 
Tstim = 5 ms have duration x0.46 x Tstim = 2.3 ms. Thus, all 
of the above cases may yield quasi-invariant spatial potentials at 
the end of stimulation, and may also be otherwise similar. 

For two representative linear-growth cases Fig. 1 1 illustrates the 
corresponding waveforms u*(t) and their construction in detail. 

Finally, Fig. 12 uses the same-vertical-scale to compare the 
relative contributions of the growth rate and the compensated re- 
polarizing node currents for each different duration. The 
waveforms' offsets (due to k*) are inversely proportional to 
duration. This readily compares qualitatively with the results in 
[8]. Especially for very short durations (e.g. Tstim =20fis), the 
optimal waveform u*{t) has a significant rectangular component 
(see also the optimality-analysis for the simple 0D models). Further 



parallels may be made for the relatively shorter durations 
(<200^). 

Numerous essential differences in the approach preclude further 
objective comparisons. Interestingly however, for the longer 
durations (T> 0.5 ms) the results in [8] show very little (if any) 
variation with Tstim (there called pulse-width, PW). 

Finally, with long PW's in [8] most of the stimulation's energy is 
delivered toward the middle of the active period. This late and 
peaky delivery requires additional analysis and comparisons of the 
actually achieved waveform-energy levels, which cannot be done 
in its details at this time. However, we return to the late delivery 
policy in the Discussion (see below), where it is deemed equivalent 
to a shorter-duration case. 

The latter provides a clue why such significant delivery 
differences would not be at odds with the very narrow 95% 
confidence intervals that resulted from the genetic algorithm in 
[8], and seeming to preclude different optimal waveforms. 

Discussion and Conclusions 

In eqn. (23), we addressed directly the electric power required for 
driving the excitable-tissue membrane potential V{t) from its 
resting (V r ) to its threshold value (Vthr) through a stimulation of 
fixed duration. Through the LAP perspective, we obtained eqn. 
(34) - a general (model-independent) description of the energy- 
optimal time-course of the excitable-tissue's membrane potential 
V(t). 

We would like to bring the reader's attention to three specific 
conclusions. 

The first is related to the intuition gained with respect to the 
evolution of the membrane potential V*(t). This optimality 
principle is best demonstrated by the simplest linear sub-threshold 
model (LM). Let ES circumstances be characterized by large 
opposing currents (e.g. the leak LM current) over long durations. 
This situation is physically analogous to filling with water a bucket 
which has large holes in its bottom. Since only the final outcome is 
important (i.e. we want the bucket full at the final time T), the best 
policy is to do nothing for most of the duration and then be able to 
dump a very large amount of water in the bucket over very short 
time. From experience, we know that works for even an unplugged 
sink. Moreover, we saw that the same intuition transfers to more 
refined models (e.g. the HHM or the MRG'02) as do nothing for 
most of the duration means that we are still around the resting V 
and hence there is no danger of Na + ionic-channel deactivation. 

The second take-home message is that the use of LAP principles 
jointly with numerical approaches (e.g. the classical FHOC) 
provides a mathematically sound and practical waveform optimi- 
zation approach, providing more assurance toward the quality of 
the final outcome. 

And finally, a note of humility is in perfect order. In this work 
we just slightly opened the door to using the LAP ideas for optimal 
ES. There are many more aspects to tackle than the ones that we 
can address in this short paper as 'proof of concept'. In particular 
we would like to extend the method for extracellular stimulation in 
forthcoming work. The motivation for doing so is at least twofold. 
On the one hand, extracellular stimulation has far more practical 
relevance. On the other hand, the only way we could rigorously 
employ the general LAP solution of eqn. (45) is to consider a 
model where we are free to manipulate u(t,Z) in every 
compartment or at every spatial location. 

A direction for such manipulation is provided by the activating 
function concept [15,20,25], which supplies every compartment 
with a virtual injected current. In the context of extracellular 
stimulation, we will also have to properly address the conditions 
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for stable AP propagation (see [15,35] for an extensive treatment 
of the subject). The optimal pattern of extracellular potentials (size 
of depolarized and hyperpolarized regions) depends on the 
distance to the electrode. These conditions would also naturally 
provide the spatial voltage profile at the end of the stimulation, 
needed to properly solve the PDE of eqn. (45). 

Here we took a shortcut path by assuming that intuitions gained 
with single-compartment models suffice. This may be partially true 
with the specific MRG'02 setup that we addressed, but does not 
hold in general. Hence, the LAP results are approximate. A clue is 
provided by the slightly lower P values of the optimal rectangular 
waveform, for TgriM =100 and 200 [is - see Table 9. As can be 
seen from Fig. 9, no benefit in terms of lower Vthr can be 
associated to the steep rise of the rectangular waveform, since 
Vthr is expected to be higher, esp. for dramatically shorter 
durations. This was further confirmed by numerical testing with 
dual linear (high/low rate) V(t) rise schedules (data not shown), 
which all had inferior performance to the baseline simple linear- 
growth protocol. However, the rectangular waveform also leads to 
steep capacitive decay of V(t) at the end of the stimulation, which 
may trigger specific patterns of additional depolarizing currents. 

For the shortest durations, the plain rectangular waveform 
outperforms by P the ones associated to the linear-ramp voltage 
profile (see Fig. 10). On Fig. 13 one can see that the steep rise of 
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the SQR waveform yields an early super-linear ramping of the 
membrane voltage. However, the rectangular waveform requires a 
lot more charge Q to be transferred. 

In practical situations many more additional aspects need to be 
addressed. E.g. stimulation needs to be charge balanced. This is a 
necessity for implanted devices and also debatably important for 
transcutaneous applications. Such stimulation will have an effect 
on the optimal threshold intensity of the cathodic pulse [36]. One 
would expect that a pre- or post- anodic pulse would also have a 
significant effect on the optimal waveform. Moreover, its own 
shape would be subject to optimization - e.g. to minimize the 
overall energy level required - a cost suitable for the design of 
implanted devices. 

We hope that the analysis and numerical evidence provided in 
this work may convince the reader of the practical benefits of 
applying the LAP principles toward the design of energy-efficient 
ES. 
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